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u,v,w  =  veloctiy  components 
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o—T* An  analytical  study  is  presented  regarding  the 
determination  of  a  pressure  gradient  dependent  criterion  for 
flow  transition  from  laminar  to  turbulent  flow.  The  results 
obtained  were  derived  from  two  parallel  approaches  to  '  ow 
stability;  one  of  wave-dependent  stability  and  the  oth  of 
vorticity-dependent  stability.  In  both  cases,  one  a 
variable  transition  Reynolds  number  dependent  upor.  ue 
ambient  pressure  gradient  and  the  other  one  of  a  constant 
transition  Reynolds  number  based  on  the  boundary-layer 
displacement  thickness,  the  prediction  results  were  either 
as  good  or  better  than  those  from  available  prediction 
methods.  In  addition  these  two  criteria  were  used  to 
predict  transition  locations  on  a  NACA  0018  airfoil,  again 


with  favorable  results.  J—yr 


v.  ,  *  «  *  .  * 
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I. 


Discussion 

The  development  of  an  understanding  of  the  mechanism  or 
mechanisms  responsible  for  the  transition  from  laminar  to 
turbulent  flow  is  one  of  the  most  unknown  areas  in  the  field 
of  fluid  mechanics  today.  Rather  than  being  dependent  upon 
one  flow  condition  or  parameter  the  actions  of  flow 
transition  is  seen  to  be  a  function  of  many  different  flow 
factors  that  are  not  necessarily  independent  of  one  another. 
Some  of  these  are:  the  roughness  of  the  flow  passage  walls; 
heating  or  cooling  of  the  passage  walls;  obs_ructions  in  the 
flow  field;  blowing  or  suction  at  the  passage  walls;  mass 
transfer  or  chemical  reaction  with  the  passage  walls; 
pressure  pulses  (noise)  in  the  flow;  deformable  wall 
geometry;  and,  ambient  pressure  gradients  in  the  flow 
(Ref.  1).  One  of  the  major  problems  in  either  an 
experimental  or  analytic  treatment  of  this  subject  is  the 
difficulty  in  isolating  the  effect  of  the  variation  of  one 
flow  parameter  while  keeping  the  others  constant.  This 
problem  is  due  to  both  the  technical  difficulties 
encountered  in  monitoring  and  controlling  the  other 
parameters  and  to  the  functional  interdependence  of  these 
parameters  which  results  in  large  data  scatters  in  similar 
investigations  and  inhibits  comprehensive  studies  of  flow 
transition.  These  complications  have  resulted  in 


significant  obstacles  in  the  theoretical  study  of  flow 
transition  and  have  limited  the  experimental  data  to  a 
small,  specialized  body.  The  many  parameters  used  to 
characterize  the  flow  take  a  wide  variety  of  functional 
forms  which  also  tend  to  defy  correlation  into  a 
comprehensive  picture. 

An  additional  barrier  to  fluid  transition  study  is  that 
not  only  are  the  factors  affecting  transition  variable  in 
their  effect  upon  transition  and  related  to  each  other  in  a 
complicated  manner,  but  transition,  and  subsequently 
transition  data,  has  not  been  uniquely  defined.  Although 
there  exist  universally  accepted  concepts  for  laminar  flow 
and  turbulent  flow,  the  transition  of  a  flow  has  been 
variously  reported  in  terms  of:  the  start  of  transition 
(the  appearance  of  turbulent  spots);  the  sud  of  transition 
(the  flow  consisting  of  a  variably-defined  level  of 
turbulence);  or,  a  region  of  transition  that  lies  between 
fully-laminar  flow  and  fully-turbulent  flow  and  contains  the 


two  previously  mentioned  boundaries.  Finally,  most 
theoretical  studies  of  this  field,  both  analytic  and 


numerical,  do  not  deal  with 


at  all  but  with  flow 


where  an  unstable  flow  is  defined  as  a  flow 
which  contains  Toll m ei n- Schl i ch t i ng  waves  of  such 
frequencies  and  wavelengths  that  their  amplitudes  will  tend 


to  magnify  as  a  function  of  time  and  will  subsequently  cause 
the  flow  transition  to  turbulence. 


Er&falsm  SisisiDfini 


With  the  above  limitations  taken  into  consideration, 
this  study  focused  on  only  the  effect  of  an  ambient  pressure 
gradient  in  the  flow  on  the  mechanism  of  transition.  An 
empirical  prediction  method  based  on  a  transition  Reynolds 
number,  based  on  fluid  boundary-layer  momentum  thickness, 
was  developed  by  Kays  and  Crawford  (Ref.  2)  that  yielded  a 
single  value  of  360  for  a  flow  regardless  of  the  history  of 
the  flow  boundary  layer.  Hermann  Schlichting  (Ref.  1), 
conversely,  reported  a  theoretical  dependence  of  fluid 
instability  on  the  presence  of  a  pressure  gradient  where  he 
used  a  critical  Reynolds  number  based  on  the  fluid  boundary- 
layer  displacement  thickness.  This  study  sought  to  develop 
a  prediction  method  that  blended  these  two  ideas  into  a 
single  criterion  (i.e.,  encompass  both  the  momentum 
interaction  of  the  boundary  layer  and  the  influence  of  the 
presence  of  a  flow  pressure  gradient). 


II.  gackground 


Sisisiilix  Uafiflci 

The  origins  of  flow  stability  theory  are  derived  from 
an  attempt  by  fluid  mechanicists  to  explain  flow  transitions 
from  first  principles.  This  theory  is  based  directly  upon 
the  incompressible,  three-dimensional  Navier-Stokes 

equations  with  the  assumption  of  one-dimensional 

u 

perturbation  velocities.  With  the  use  of  pert  rbation 
theory,  the  assumption  of  small  amplitude  disturbances,  and 
the  retention  of  only  the  linear  perturbation  terras,  the 
Navier-Stokes  equations  are  transformed  into  what  is 
commonly  known  as  the  Orr-Sommerf eld  equation.  A  full 
derivation  of  the  Orr-Sommerf eld  equation,  based  in  part 
upon  the  treatment  of  Hunt  (Ref.  3)  and  in  part  upon  that  of 
Schlichting  (Ref.  1)  is  given  in  Appendix  A.  The  usual  form 
of  the  equation  is: 

( u-c) ( 0"  -  o20)  -  u"0  =  —  ( 0”  -  2<*20"  +  “40)  (1) 

aR 


where 


u  =  boundary  layer  velocity  profile 
c  r  complex  wave  speed 
0  s  wave  amplitude 


a  =  wave  number 


f  / 


(J'J'Ji  TOR  'I'.'JL  P  H  W  'A  A  '/• 7 


i  =  yfT 

R  =  Reynolds  number  (usually  based  upon  edge 

velocity  and  displacement  thickness) 


A*  > 

V. 

J"  • 


and  all  variables  are  nondimensional  with  the  primes 
denoting  derivatives  with  respect  to  the  coordinate  normal 
to  the  mean  flow  direction. 

The  Orr-Sommerf eld  equation  represents  a  complex, 
fourth-order  homogeneous  ordinary  differential  equation  in 
0,  yet  with  both  u  and  0  being  functions  of  the  normal 
coordinate  and  due  to  the  fact  that  (C, a  ,  R)  represent  a 
related  family  of  parameters,  there  is  no  closed-form 
solution  to  this  equation.  Thus,  the  main  tool  of  flow 
stability  has  no  known  exact  solution.  Common  approximate 
analytic  solutions  use  asymptotic  approximations  that 
involve  recursive  Hankel  functions  of  order  1/3,  such  as 
Lin's  method  (Ref.  4).  Existing  computer-aided  numerical 
solutions  as  a  rule  either  use  eigen-value  "shooting 
methods"  (Ref  5)  or  expand  the  equation  in  terms  of 
orthogonal  polynomials,  such  as  Chebyshev  polynomials 
(Ref.  6),  and  use  su c c e s s i v e - gu e s s  and  recursive 
relationships  among  the  polynomials  to  achieve  approximate 
solutions  to  the  Orr-Somraerf ield  equation. 

The  usual  method  of  presenting  these  stability  results 
is  in  the  form  of  a  neutral  stability  curve  (see  Fig.  1). 
In  this  plot,  the  curve  represents  the  variation  of  wave 
number  as  a  function  of  Reynolds  number.  The  points  on  the 
curve  represent  Reynolds  numbers  and  wavenumber  combinations 
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Pig  1  A  Typical  Stability  Diagram 


which  do  not  tend  to  be  amplified  or  dampened  temporally  for 


the  associated  waves.  The  region  outside  of  the  curve 
represents  those  combinations  of  wavenumber  and  Reynolds 
number  for  which  disturbances  tend  to  damp-out  and  the 
region  inside  the  curve  represents  those  combinations  for 
which  disturbances  tend  to  be  amplified  and  for  which  the 
flow  is  said  to  be  unstable.  The  Reynolds  number  below 
which  no  value  of  wavenumber  will  yield  unstable  flow  is 


said  to  be  the 


Reynolds  number  and  is  the  minimum 


value  of  the  Reynolds  number  on  the  curve.  Due  to  the 
approximate  nature  of  the  existing  solutions  there  is  a 
considerable  variation  of  determined  Rcrit  values  as  can  be 
seen  in  Fig.  2  which  is  from  a  compilation  by  van  Ingen 
(Ref.  7)  where  Rcrit  varies  from  321  to  680.  With  this 
large  variation  in  mind,  stability  theory  results  are 
treated  warily  and  are  considered  as  approximate  only. 


a 


Stability  theory,  though  it  is  derived  from  first 
principles  through  the  Navier-Stokes  equations,  is  in  fact 
quite  limited  when  true  transition  is  considered.  Besides 
the  approximate  nature  of  the  Orr-Sommerf eld  equation 
solution  which  does  not  yield  unique  values  for  any  one  flow 
situation  treated  by  stability  theory,  most  common  flows 
tend  to  violate  one  or  more  of  the  basic  assumptions  used  to 
derive  the  theory.  Some  of  these  violations  are:  velocity 
components  nonparallel  to  the  mean  flow  direction  that  tend 


Elg  2  Neutral  Stability  Curves  for  Flat  Plate  Flow 
with  Zero  Pressure  Gradient  from  Different 
Sources.  (Ref. 7) 


to  make  the  flow  unstable  at  lower  Reynolds  numbers; 
nonlinearity  of  the  perturbation  velocity  components  that 
tend  to  invalidate  the  linearized  theory  used  to  derive  the 
Orr-Sommerf eld  equation;  dependence  of  the  mean  flow 
velocity  on  its  streamwise  coordinate  that  is  neglected  in 
stability  theory  derivations  and  adds  another  degree  of 
complexity  to  the  flow  situation;  and,  large  amplitude 
disturbances  that  are  observed  just  prior  to  transition  that 
are  considered  infinitesmal  in  magnitude  for  stability 
theory.  Even  with  these  severe  limitations  of  stability 
theory,  Rcrit  is  always  less  than  Rtrans  for  a  given  flow 
case  and  thus  provides  a  conservative  criterion  for 
transition.  Furthermore,  the  behavior  of  Rcrit  from  one 
flow  case  to  another  tends  to  follow  the  same  trends  that 
Rtrans  does,  yielding  a  characteristic,  if  not  exact, 
criterion  for  flow  transition. 

The  fundamental  basis  of  flow  stability  theory  is  that 
flow  instability,  and  thus  flow  transition,  is  dependent 
upon  the  behavior  of  the  amplification  or  dampening  of  the 
Tollmein-Schlichting  waves  present  in  the  flow.  A  mechanism 
that  does  not  depend  on  this  type  of  wave  action  yet  could 
plausibly  play  an  equal  if  not  greater  role  in  flow 
transition  would  be  a  welcome  addition  to  transition  theory 
if  it  could  be  seen  to  have  a  firm  foundation  in  first 
principles  and  also  not  to  have  many  of  the  limitations  of 
wave-theory.  Such  an  approach  which  was  based  on  "vorticity 
stability"  was  investigated  in  this  study  with  favorable 
results. 
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III.  Theory  Development 


The  determination  of  a  pressure  gradient  criterion  for 
flow  stability,  and  by  allusion  transition,  required  the 
knowledge  of  the  values  of  A(Pohlhausen  velocity  profile 
shape  parameter),  6  (boundary  layer  thickness),  51 
(displacement  thickness) , and  5  2  (momentum  thickness)  for 
each  flow  case  in  addition  to  the  value  of  Rcrit.  This  was 
due  to  the  fact  that  various  sources  based  their  cited 
values  of  Rcrit  or  Rtrans  on  either  51  or  «2,  or  in  the 
case  of  internal  flow,  Dh  (hydraulic  diameter)  which  could 
be  related  physically  to  for  fully-developed  internal 
flows. 


Emaaurfi  Gradient 


The  pressure  gradient  parameter  that  was  used 
throughout  this  study  was  Pohlhausen's  shape  factor  based  on 
boundary  layer  thickness, A  ,  or  a  modified  form  commonly 
used  which  is  based  on  momentum  thickness,  x  .  The  normal 
Pohlhausen  parameter  is  defined  as 


A 


5^  dUm 
v  dx 


(2) 


where  Urn  =  Um(x),  and  in  the  case  of  the  wedge  family,  Urn 
is  of  the  form,  Urn  =  Uoxm,  with  Uo  being  a  constant  and  x  is 
the  streamwise  coordinate,  x  is  defined  similarly  as 


x  _  2  dUm 


For  the  flow  cases  which  were  considered  here,  the  flow  was 
attached  laminar  flow,  i.e.  instability  and  transition 
occured  before  flow  separation  from  the  boundary  wall.  In 
the  case  of  external  flow,  Bernoulli's  equation  for  along  a 
streamline  outside  the  boundary  layer  is  given  as 


dP  _  dUm 

dx  dx 


and  this  can  be  used  with  equation  (1)  to  denote  a  and  \  as 
pressure  gradient  parameters  of  the  form: 


,  (i!) 

yUm  ' dx  ' 


X  =  -mil  /") 

yUm  ' dx7 


This  rendering  of  as  a  measure  of  the  magnitude  and 
typeof  pressure  gradient  (  a  positive  denoting  a  streamwise 
pressure  drop,  or  favorable  pressure  gradient,  etc.)  was 
very  convenient  for  the  flat  plate  cases  studied  due  to  the 
common  treatment  (Ref.  1)  of  the  boundary  layer  profile  as  a 
polynomial  in  A .  Here,  the  velocity  profile  was  of  the 
form: 


*  m  *  ^  «  »  m  *_»  "  m  “  —  m  *  jiT  m  ^  »  m  j, 


(7) 


u(  n)  s  Urn  - ( 2n-2  n3+ )  +  —  ( ri-3n  2+3n3- ) 

6 

where  n is  a  nondimensional  coordinate  to  the  flow  commonly 
denoted  as  ns  y/«.  With  the  variation  of  A  t  different 
profiles  of  the  nondimensional  (corresponding  to  u(n)/Um) 
in  equation  (1))  were  obtained  by  Schlichting  (Ref.  1),  et. 
al.,  that  yielded  different  Rcrit  values  depending  upon  the 
value  of  a  (see  Fig.  3).  This  functional  relationship  of 
Rcrit,  based  on  displacement  thickness,  was  plotted  versus  a 
in  Fig.  4  (the  plot  as  far  as  can  be  determined  from  Ref.  1 
and  Ref.  8  was  merely  a  curve  drawn  through  the  indicated 
stability  trials).  Schlichting  further  implied  that  this 
relationship  for  Rcrit  vs.  a  ,  based  on  various  pressure 
gradients  imposed  on  flat  plate  flow,  could  be  universally 
applied  to  all  external  boundary  layer  flows. 

The  assumption  of  a  universal  relationship  between 
Rcrit  and  a  as  given  in  Fig.  2  was  tested  by  considering 
the  "wedge  flow"  cases,  which  are  particular  solutions  of 
similar  solutions  of  boundary  layer  flow.  These  flows 
experience  a  streamwise  pressure  gradient  that  is  due  to  a 
flow  area  change  where  the  conventional  parameter  is 
Hartree’s  8 .  Hartree’s  8  is  physically  interpretted  as  the 
ratio  of  the  wedge  angle  over  and  is  a  function  of  the 
exponent  of  the  stream-wise  coordinate  for  similar  flows 
denoted  in  the  following  form: 
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Fig  3  Neutral  Stability  Curves  for  Flat  Plate 
with  Pressure  Gradients.  (Ref.  1) 


z  ♦  t  t 


Fig  4  Critical  Reynolds  Number  Variation  with 
Pressure  Gradient  Factor,  A.  (Ref.  O 


m  +•  1 


The  relationship  between  a  and  e  (see  Appendix  B  for  a 
full  derivation)  was  found  to  be 

a  =  — (  n.)2  B(m  +  1 )  (9) 

2  6 

where  n6  is  the  value  of  the  nondimensional  coordinate  at 
the  edge  of  the  boundary  layer  (where  u(n)  is  99%  of  Um). 
With  this  relationship  it  was  possible  to  compare  on  one 
plot  (see  Fig.  5)  both  the  flat  plate  cases  and  the  wedge 
flow  results  and  to  verify  that  what  Schlichting  implied  but 
did  not  prove  appeared  to  be  correct. 

The  two  internal  flow  cases  of  two  parallel  infinite 
flat  plates  (plane  Poiseuille  flow)  and  fully-developed  pipe 
flow  (Hagen-Poiseuille  flow)  were  also  considered  in  this 
study.  There  was  a  stability  result  available  for  the 
parallel  flat  plate  case  (Ref.  6)  but  the  pipe  flow  case  was 
found  to  be  stable  to  inf intesiraal  disturbances  up  to  all 
critical  Reynolds  numbers  tried  (Ref.  9)  and  no  stability 
result  was  thus  available.  However,  a  corresponding 
approximate  parameter  was  estimated  for  both  cases  due  to 
the  availability  of  transition  data  for  both  cases.  This 
approximate  method  followed  from  the  known  exact  solutions 
of  the  Navier-Stokes  equations  (Ref  1:85-86)  for  Poiseuille 
flow  of 

-h2  dP  y  2 

u(y)  =  -  (—Ml  -  r-r)  (10) 

2v  dx  n 
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and  for  Hagen-Poiseuille  flow  in  a  pipe  given  by 

-R2  dP  y  2 

u ( y )  =  -  ( — ) ( 1  -  (-)  )  (11) 

4  u  dx  R 

where  y  is  a  coordinate  measured  from  the  centerline,  x  is 
the  streamwise  coordinate,  h  is  half  the  distance  between 
the  parallel  plates,  and  R  is  the  radius  of  the  pipe.  The 
relationship  between  the  pressure  gradient  and  from 
equation  (5),  with  the  substitution  of  U  centerline  for  Um 
and  h  and  R  respectively  for  5,  was  adapted  wholesale  to 
obtain  "equivalent”  values  of  A  for  the  parallel  plates  and 
the  fully-developed  pipe  flow,  which  were  respectively, 
a  =  2  and  a  =  4.  The  values  of  a  obtained  were  those 
necessary  to  be  consistant  with  the  parabolic  velocity 
distributions  commonly  given  for  these  flow  cases. 

Ifalgfcngaasa 

For  the  external  flow  cases  the  different  boundary- 
layer  thicknesses  were  defined  in  the  usual  manner.  The 
usual  boundary-layer  thickness,  s  ,  was  the  distance  out  from 
the  wall  at  which  the  local  velocity  had  obtained  99%  of  the 
"free- stream "  velocity  value.  The  displacement  thickness, 
1 ,  was  given  by 

f  U(n) 

«1  =  /  (1 - )  dn  (12) 

o  U„ 

And,  the  momentum  thickness,  62,  was  given  by: 


f  U(n)  ,,  U(n1 

5  2=  /  -  ( 1  - - )d  n  (13) 

o  U»  U« 

The  flow  thicknesses  for  the  interior  flow  cases  were 
obtained  by  modifying  the  conventional  definitions  given 
above.  The  boundary  layer  thickness,  6  ,  was  assumed  to  be 
half  the  flow  passage  dimension  for  both  the  fully-developed 
internal  flow  cases.  A  detailed  analysis  of  the 
approximations  for  51  and  52  of  both  internal  flow  cases  is 
given  in  Appendix  C.  The  values  for  the  parallel  flat 
plates  were  approximately  51  =  h/3  and  s2  =  h/7.5  where  h 
was  half  the  plate  separation,  and  the  values  for  the  fully- 
developed  pipe  flow  were  approximately  51  s  (0.293)R  and 
52  =  (0.1 835 ) R »  where  R  was  the  pipe  radius. 


The  determination  of  the  approximate  value  of  the 
displacement  thickness  for  the  parallel  plates  case  along 
with  its  "equivalent"  Pohlhausen  parameter  allowed  the 
inclusion  of  another  stability  theory  result  (Ref.  6)  on  the 
plot  representing  the  variation  of  Rcrit  with  a (see  Fig.  5). 
As  can  be  seen  this  plot,  the  internal  flow  result  seem  to 
be  consistent  with  the  external  flow  cases  for  this  work's 
interpretation  of  a  . 

Since  most  transition  data  appeared  to  be  reported  in 
terms  of  R52  and  not  R  5 1 ,  the  stability  results  were 
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converted  into  critical  Reynolds  numbers  based  on  «2  (see 
Fig.  6).  An  analytic  equation  was  then  determined  that  fit 
the  stability  results  that  was  given  by 


A  -  4.5 

R«?»crit  =  2954  tanh  ( - — )  +  2981 

2.7 


(14) 


and  was  presented  with  the  stability  results  in  Fig.  7. 
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WEDGE  FLOWS 
PARALLEL  PLATE  CASE 
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Fig  8  Comblnad  Stability  Raaulta. 
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IV.  SJahll Hi 


Concurrently  with  the  work  done  to  produce  an  analytic 
expression  for  the  wave  stability  dependence  on  the  ambient 
pressure  gradient,  the  stability  of  the  vorticity  of  the 
flow  was  also  investigated.  The  laminar  flows  considered  in 
this  work  were  assumed  to  be  parallel.  General  turbulent 
flow,  however,  is  known  to  be  greatly  unparallel  locally 
where  large  vortices  feeding  "large  eddies"  break  up  the 
main  flow  and  in  turn  devolve  into  "small  eddies"  that  are 
responsible  for  turbulence  on  a  local  level  (Ref.  10).  As 
long  as  the  viscous  stresses  acting  against  the  flow 
counter-balance  the  flow’s  inclination  to  rotate  and  form 
eddies,  the  main  mechanism  of  turbulent  flow  motion,  eddy 
production,  cannot  be  sustained  and  the  flow  remains 
laminar.  The  objective  of  the  work  done  in  this  section  was 
to  determine  a  means  of  predicting  when  the  flow  would 
"trip"  over  into  eddy  motion  (rotating  pockets  of  fluid)  and 
thus  become  turbulent. 

The  anlaysis  of  the  effect  of  vorticity  upon  flow 
transition  was  based  upon  a  rough  visualization  of  the  flow 
as  a  solid  element  moving  along  a  resisting  surface 
(Fig.  8).  As  long  as  the  moment  due  to  the  element’s  weight 
counters  the  inertial  moment,  the  element  will  tend  only  to 
translate  along  the  contact  surface.  If  however  the  element 
is  elongated  away  from  the  contact  surface,  the  weight 
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moment  will  decrease  and  the  inertial  moment  will  increase 
(Fig.  8b).  If  this  deformation  is  continued  until  the 
inertial  moment  is  greater  than  the  weight  moment  (Fig.  8c), 
the  element  will  tend  to  rotate  as  well  as  translate.  When 
this  solid  element  is  now  considered  as  an  incompressible 
fluid  element,  the  inertial  moment  corresponds  to  the 
product  of  the  shear  stress  (proportional  to  the  vorticity) 
in  the  fluid  and  its  location  away  from  the  flow  passage 
wall.  The  flow  (which  is  the  summation  of  all  of  these 
infinitesimal  moving  fluid  elements)  is  considered  to  be 
unstable  when  the  first  moment  of  the  flow  vorticity  is 
greater  than  the  effect  of  the  net  balancing  forces  acting 
on  the  fluid  element  (Ref.  11).  While  these  balancing 
forces  are  not  identified  specifically,  it  is  hypothesized 
that  they  are  dependent  on  fluid  properties  but  independent 
of  local  flow  conditions;  that  is  to  say,  once  the  first 
moment  of  vorticity  gets  to  be  some  magnitude,  the  flow 
becomes  unstable. 

The  vorticity  of  the  parallel  flows  considered  was 
given  by 

c  =  V  X  V  =  (—  -  — )  (14) 

ax  ay 

which,  if  av/ax  is  neglected,  was  simplified  to 

c  =  -  —  (15) 

ay 

The  first  moment  of  vorticity  was  designated  as  I,  where  I 
was  defined  by 
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=  -  /  -  cydy 


(16) 


The  viability  of  this  approach  was  tested  by 
considering  the  Pohlhausen  boundary-layer  profile  and  the 
parabolic  velocity  profile  for  fully-developed  parallel 
plate  flow.  In  the  first  case,  the  velocity  profile  was 
given  by 


u(y)  .  u  -  2z!  ♦  z!,*  4  (L.  -  -  £->>  (17) 


6  5£ 


fi3  «4 


The  moment,  I,  was  then  given  by 


O 

U  C 
«  -  /  < 


2y  6y3  4ya 


«  «■ 


A (y  6y2  9y3 

6  a  “  *2  +  53 


4y4 

— jp))  dy 


with  the  result  that 


i  =  Hi(J  .  _jl) 

V  MO  1 2  0 y 


(19) 


Using  the  shape-factor  ratio  of  ^  to  6  given  by  Pohlhausen 
(Ref.  1)  as  (3/10  -  a/120),  this  measure  of  stability  was 
then  I  s  R  5  .j  f  or  the  Reynolds  number  basedon  the 
displacement  thickness.  The  second  case  used  a  velocity 
profile  given  by 
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where  half  the  displacement  between  the  plates,  h,  was 


equated  to  the  boundary-layer  thickness,  6  .  The  moment,  I, 
was  then  given  by 

6 

U  r  2y  2y2 

I  =  _  /  (—  -  JL)  dy  (21) 

v  J  & 

O 

and  resulted  in, 


When  compared  with  the  shape-factor  previously  determined 
for  the  parallel  plates  case  of  1/3,  I  again  was  the  value 
of  the  Reynolds  number  based  on  the  boundary-layer 
displacement  thickness,  R«-j  .  These  results  indicated  that 
the  value  of  the  Reynolds  number  based  on  displacement 
thickness  was  the  single  parameter  for  flow  stability. 


The  validity  of  the  two  prediction  methods  was  examined 
by  comparing  the  predicted  results  with  those  available  from 
experimental  investigations  in  two  different  ways.  The 
prediction  results  were  first  compared  to  accumulated  data 
for  a  variety  of  different  airfoils  and  flow  cases  and  then 
the  prediction  methods  were  used  to  analytically  estimate 
the  location  of  the  transition  point  for  a  NACA  0018  airfoil 
that  experienced  two  different  body  Reynolds  numbers  and 
different  angles  of  attack. 

frsaaurg  SradlgD-fc  Method 

The  pressure  gradient  method  that  was  derived  from 
stability  theory  was  converted  into  transition  data  in  two 
different  ways.  Initially  the  analytic  curve  that  fit  the 
stability  data, 

R«2,crit  5  2954  tanh  +  2981  (23) 

was  converted  to  a  transition  basis  by  scaling  the  predicted 
critical  value  of  852*crit  =  231  for  the  zero  pressure 
gradient  case  (As  0)  to  the  transition  value  of 

R«2*trans  =  1150.  This  transition  value  was  obtained  by 
modifying  the  Kays  and  Crawford  method  (Ref.  1)  where  a 
value  for  Rx,trans  =  3,000,000  was  used  instead  of  their 
value,  Rx,trans  =  300,000  (this  higher  value  is 
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experimentally  supported  by  Reference  12).  The  resulting 
analytic  equation  was  then: 


R«2»trans  =  14786  tanh  +  14917  (24) 


The  second  method  took  into  account  the  known  difference 

between  Rs2,crit  and  R«2»trans  shown  in  Figure  9  (Ref.  13). 
An  analytic  curve  fit  for  this  difference  was  found  to  be 


aR$2  =  400  exp(60 x)  +  400 


(25) 


where  x  was  related  to  a  (Ref  2:210)  by 


x  ,  ( Jl .  _l_A _ L 

315  945  9072 


a2)2, 


(26) 


This  difference  was  then  added  to  the  R«2»crit  expression  of 
equation  (23)  to  obtain  a  transition  Reynolds  number  of  the 


form : 


R$2>  trans  =  2954  tanh  )  +  2981  +  ar52  (27) 


The  results  of  both  of  these  methods  are  shown  in  Figure  10, 
together  with  experimental  results  which  were  taken  from 

Reference  14.  (It  must  be  noted  that  the  expression  for  AR $2 
was  only  valid  for  x  values  less  than  0.025,  or  a  values 
less  than  2.0,  and  subsequently  care  should  be  taken  when 
using  that  method  of  conversion.) 


♦  •  EXPERIMENTAL  TRANSITION  DATA 
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Fig  10  Comparison'  of  Stability— 
Derived  Curves  with  Data. 
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l&Lt. Lziti  SiaMliii  Heiiaad 


Two  values  of  Rfiiitrans  =  constant  were  chosen 
to  represent  the  vorticity-stability  theory;  a  value  or 

R«1  =  3000  that  corresponded  to  low  (0%)  free-stream 
turbulence  levels  and  a  value  of  R  ^  =  1226  that 

corresponded  to  high  ( 1 .2 % )  free-stream  turbulence  levels. 
These  values  were  chosen  so  as  to  be  consistent  with  the 
findings  of  Hall  and  Gibbings  (Ref.  13).  A  shape-factor 
consistent  with  Pohlhausen’s  polynomial  fit  (Ref.  1:210)  was 
used  to  relate  R$  ^  to  Rs-j  by  the  function 


R52  =  R61 


37 _ ]_ _ 1_ 

315  945  A  9075  A 


(28) 
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The  resulting  variations  of  R«2>trans  with  x  are  shown  in 
Figure  11  along  with  previous  results  and  further  data  (from 

Ref.  13).  It  was  evident  that  a  constant  value  of  Rg  i  fit 
the  data  for  x  <  0  (adverse  pressure  gradients)  quite  well 

but  tended  to  underestimate  the  value  for  Rs2»trans  **or 
x>  0  (favorable  pressure  gradients).  This  lack  of 
correlation  for  vorticity  stability  to  experimental  data  for 
accelerated  flow  might  be  due  to  the  limitations  of  a  fourth 
order  polynomial  fit  for  the  velocity  profile  in  this  region 
or  might  be  due  to  other  undetermined  reasons.  It  is 
interesting  to  note,  however,  that  in  the  favorable  pressure 
gradient  region  the  transition  equations  based  on 
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wave-stability  fit  the  low-level  free-stream  turbulence 
data.  It  is  possible  that  wave  mechanics  effects  over¬ 
shadow  the  effects  of  vor t i ci ty- sta bil i ty  of  the  flow  for 
these  cases. 

The  two  values  of  R{1  =  constant,  3000  and  1226,  that 
are  presented  in  Fig.  11  are  arbitrary  bracketing  cases  for 
the  data.  These  values  in  turn  represent  flat  plate 
Reynolds  numbers  based  on  the  location  from  the  leading  edge 
of  363,000  and  3,130,000  respectively,  which  are  within 
reported  experimental  limits  of  experimental  flat  plate, 
zero  pressure  gradient,  transition  (Ref.  12). 

Internal  Elau  flaaalte 


The  two  internal  flow  results  considered  in  this  work 
concerned  parallel  flow  between  two  infinite  plates  and 
fully-developed  flow  in  a  pipe.  The  parallel  plates  result 
was  based  on  a  thin  rectangular  duct  study  where  the  aspect 
ratio  was  8:1  (width/depth)  and  the  transition  Reynolds 
number  based  on  average  velocity  and  hydraulic  diameter  was 
experimentally  found  to  be  2600  (Ref.  15).  This  Reynolds 

number  yielded  R$2*trans  =  58s*  The  pipe  flow  result  was 
based  on  the  universally  quoted  value  of  the  Reynolds 
number,  based  on  hydraulic  diameter  and  average  velocity,  of 
2000  for  transition.  This  value  yielded  a  value  of 
R52,trans  s  both  of  these  results  are  shown  in  Fig.  11 
and  signified  by  "o".  It  is  interesting  to  note  that  the 
low-level  free-stream  turbulence  value  of  R $i  =  constant 
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(3000)  yielded  a  Rd(trans  =  6853  f  which  is  well  within  the 
known  experimental  range  of  2000  to  40,000  (Ref.  16). 

EcfidigMfifl  £X  Irsniiiifin  Locaiisn 

The  prediction  methods  developed  in  this  work  were  used 
to  estimate  the  location  of  transition  on  a  NACA  0018 
airfoil.  Schlichting  (Ref.  1:500)  gave  a  presentation  where 
the  experimental  transition  curves  of  two  body  Reynolds 
numbers,  those  of  1,700,000  and  5,000,000,  for  the  NACA 
airfoil  were  compared  with  stability  curves  of  a  Joukowski 
0015  airfoil  that  had  a  similar  pressure  distribution.  The 
work  presented  here  used  a  computer  program  (which  was 
obtained  from  Captain  J.  Lawrence,  Ref.  17)  that  calculated 
boundary-layer  parameters  at  incremental  steps  along  a 
Joukowski  airfoil.  At  various  angles  of  attack,  estimates 
were  then  made  to  locate  the  transition  points  on  the 
airfoil  for  the  various  angles  of  attack  and  resultant 

coefficients  of  lift.  Transition  curves  for  both  a  R^  -j 
value  of  1355  and  the  scaled  pressure  gradient  dependent 
method,  Equation  (24),  were  used  to  produce  the  results 
presented  in  Figures  12  and  13.  It  can  be  seen  in  both 
cases  that  the  predicted  transition  curves  did  a  better  job 
at  predicting  transition  than  the  stability  curves  alone. 
The  constant  R$i  curve  gave  a  closer  approximation  to  the 
experimental  curve  for  the  Rx  s  1.7  million  case  but  did  not 
do  as  well  as  the  pressure  gradient  dependent  method  for  the 
=  5.0  million  case. 
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VI.  Conclusions  and  Recommendations 


This  work  produced  two  different  methods  for  predicting 
the  transition  of  laminar  to  turbulent  flow.  The  first 
approach  resulted  in  a  pressure  gradient  dependent  technique 
that  matched  the  available  experimental  data  well  in  the 
region  of  adverse  pressure  gradients  but  underestimated  the 
data  for  favorable  pressure  gradients.  The  second  method, 
derived  from  viewing  the  stability  of  the  flow's  vorticity, 
resulted  in  a  criterion  of  a  constant  Reynolds  number,  based 
on  boundary-layer  displacement  thickness,  that  fit  the  data 
well  for  favorable  pressure  gradients  but  underestimated  the 
data  for  adverse  pressure  gradient  cases. 

Future  work  in  this  area  could  include  a  more  rigorous 
examination  of  the  modelling  techniques,  a  high-order 
polynomial  fit  for  the  velocity  profile,  e.g.,  that  might 
improve  the  prediction  results  for  both  methods.  In 
addition  Hall  and  Gibbings  (Ref.  13)  noted  that  Rsi,trans 
varied  with  free-stream  turbulence  levels.  This  variation 
could  be  investigated  further  or  the  effect  of  other 
parameters,  such  as  heat  transfer,  upon  R  § -j  could  be 
examined. 


N.  V  M  \  \  \ 


Afifisndii  A 

Orr-Sommerf eld  Equation  Derivation 

The  Orr-Sommerf eld  equation  is  based  on  the  three- 
dimensional,  nonsteady,  encompressible  Navier-Stokes  equa¬ 
tions,  as  follows: 
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Simpliflying  assumptions  are  that  body  forces  are  negligible 
and  that  there  is  only  one-dimensional  mean  flow  in  the 
x-direction.  The  total  flow  variables  are  sums  of  mean  and 
perturbation  quantities  of  the  form: 

u  =  u(y)  +  u' 
v  =  v' 
w  =  w  ’ 

P  =  P  +  P* 
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These  variables  are  substituted  back  into  the  Navier-Stokes 


equations  with  only  the  terms  linear  in  u',  v',  or  w' 

retained.  The  equations  are  now  of  the  form: 


u* t  +  UU'X  +  V • Uy 

P  X 

=  -  -  +  V  V2u  ’ 

p 

( A5) 

v  1 1  +  uv ' x 

p,y  2  , 

p 

(A6) 

W' t  +  uw *  x 

P'z  5 

=  -  -  +  V  v^w ' 

p 

(A7) 

U'x  +  V  '  y  +  W'2 

=  0 

(A8) 

The  disturbances  are  considered  to  be  bounded  as  x  and 
z  tend  toward  infinity.  They  are  also  considered  spatially 
periodic  of  the  form 

u'  =  u(y , t)exp[i(kx  +  lz)3 
V  =  v(y , t)exp[i(kx  +  lz)] 
w'  =  w(  y , t ) exp  t i ( kx  +  lz)] 

P'  =  P(y, t)exp[i(kx  +  lz)] 

The  equations  are  now 

^  ^  p 

ut  +  u(iku)  +  vUy  =  -  ik  -  +  v(-k2  -  l2  +  uyy)  (A9) 

p 

Py 

Ut  +  u(ikv)  =  -  +  v(-k2  -  l2  -  vyy)  (A10) 

wt  +  u(ikw)  - - +  v(-k2  -  12  _  w  )  (ah) 

p 

iku  +  Vv  +  ilw  =  0  (A12) 


Equation  A9  is  multiplied  by  k/(k2  +  l2)1'2  and  equation  All 
is  multiplied  by  l/(k2  +  I2)"*/2  and  the  two  are  added 


together  to  yield: 


-rCku  +  lw)  +  u(ik)(ku  +  lw)  -  uvku  = 
at  y 


i(k2  +  l2)-  +  v[-k2(ku  +  lw) 


-  l2(ku  +  lw) 


+  — x  (ku  +  lw) ] 
ay^ 


(A13) 


New  variables  are  now  defined  as 


U  =  ,  -2  s  k2  +  12  ,  P  =  *  ,  t  =  ^ 

(k2  +  l2)1/2  k  ® 


where  the  equations  are  now  of  the  form: 


+  u(iou)  +  vuy  =  -ia  ~  +  —  (-«2  +  uyy) 


(A14) 


't  +  u(iov)  =  -a-~  +  (-a2  +  vyy) 


CA15) 


iaUY  +  yv  =  0 


( A1 6 ) 


Equations  A14-A16  are  now  of  the  form  of  equations  A9-A12  if 


w  and  1  are  zero,  which  is  the  form  of  the  two-dimensional 


perturbation  case.  Since  (a/k)v  =((k2  +  l2)/k)v  >v  ,  then 


^crit-2d  <  Rcrit-3d  and  two-dimensional  perturbations  are 


less  stable  than  three-dimensional  ones.  This  result  is 


known  as  Squire's  theorem.  Further  analysis  can  be 


simplified  conservatively,  therefore,  by  considering  only 


two-dimensional  perturbations. 


rw 


If  equations  A9,  A10,  and  A12  are  taken  with  w  =  0 
the  mean  flow  equations  are  assumed  to  satisfy  the 
equations,  i. e. , 


i 

P  3  X 


d2u  ,  1  aP 

m  -  and  -  — 

dy2  p  ay 


0 


then  the  N.S.  equations  are  then  of  the  form: 

u't  +  uu'x  +  v*uy  +  1/p  P’x  =  w2u' 

V't  +  UV  1  x  +  1/P  P'y  =  \)72V 
u»x  +  v'y  =  0 

A  stream  function,  *  ,  is  now  introduced 
*(x»y»t)  =  0(y)exp[i(ax  -  st)] 

where 

0  =  complex  amplitude 
2 

a  r  -J  (wave  number) 

6  s  8r  +  iBi 

er  =  circular  frequency  of  partial  oscillation 
Sis  amplification  factor 
C  =  e/a 

Cr  -  velocity  of  wave  propagation  in  x-direction 
s  degree  of  amplification/damping 


The  perterbation  velocities  are  now  defined  by 


=  0* (y)exp[i( 3x  -  et) ] 


and 

N.S. 


(AIT) 
( A 1  8 ) 

(A19) 
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=  -ia0(y)exp[i( ax  -  st)] 


substituting  back  into  A17,  A18,  and  A19  yields 

0(-ie)exp[i( ax  -  e t ) 3  +  u0 1 iaexp[i( ax  -  st)] 

Px 

-  ia0exp[i(ax  -  et)3u'  +  — 

p 

=  v  ( -0 '  a^exp  [  i(oX  -  St)]  +  0  '  '  ’  exp[i  (  aX  -  fit)])  ( A20  ) 

_  _  Py 

-o80exp[i(ax  -  St)]  +  Ua20exp[i(aX  -  st)]  +  — 

p 

s  v(iUo30exp[i(aX  -  st)]  -  io0' ' !exp[iax  -  Bt)]  (A21) 

0 ' ( ia )exp [ i (a x  -  Bt)]  -  (io€'exp[i(ax  -  Bt)]  =  0  (A22) 

If  equation  A20  is  differentiated  with  respect  to  y, 
Equation  A21  is  differentiated  with  respect  to  x,  the 
pressure  term  is  assumed  to  be  indifferent  to  order  of 
differentiation,  the  resultant  Equation  A21  is  subtracted 
from  the  new  Equation  A20,  and  the  resultant  combined 
equation  is  divided  through  by  exp(i(ax  -  st)),  then  the 
single  combined  equation  is  of  the  form: 

i  80"  +  iau0"  -  ia0u"  +  ia2s0  -  ia3u0  = 

v(0«»"  -  20"  «2  +  a^0)  (A23) 


And  finally,  if  Equation  A23  is  divided  through  by  i a  (where 


L  =  8  / a)  ,  nondimensionalized  by  R  =  lUm/v,  and  the  terms 
on  the  left-side  of  the  equation  rearranged,  the  common  form 


of  the  Orr-Sommerf eld  equation  is  produced  as: 


(u  -  c)(0"  -  a2g)  »  u”0 


A£££Mi£  £ 

Conversion  of  8  Parameter  To  a  Parameter 


Pohlhausen's  shapefactor  is  usually  defined  as: 
a  _  *£  Mm 

v  dx 


(B1 ) 


where 

6  =  boundary-layer  thickness 
v  =  fluid  kinematic  viscosity 

Um  =  (U0  =  constant) 

x  =  streamwise  coordinate 


Hartree's 
defined  as 

e  = 


beta  parameter 
(Ref.  1): 

L  2  dUffl 
U.  •  8  dx 


or  wedge  angle  (see  Fig. 


5)  is 


(B2) 


where 

L  =  characteristic  length 
=  free-stream  velocity 

and 


g 


m  +  1 


>  'l  '5 


m 


(B3) 


If  Equation  (B3)  is  substituted  into  Equation  (B2),  the 
resultant  expression  for  Hartree's  parameter  is 


8 


m  +  1 


) 


x 

U 


dUra 

w 


( B4 ) 


’  -  v  v 
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Equations  (B1)  and  (B4)  are  then  solved  for  dUm/dx  t< 
yield: 


8 .  Um  .  (m  +  1 )  Av 
2  x  62 


(B5) 


Rearranging  yields, 


Um52 

A  =  (* - -)  .  0  .  (m  +  1 ) 

d  .  v  x 


(B6) 


For  a  given  boundary-layer  flow,  the  boundary-layer 
thickness  is  given  in  terms  of  the  similarity  variable,  n  , 
in  the  form: 


6  =  n, 


where  n5  is  the  value  of  the  similarity  variable  at  the  edge 
of  the  boundary  layer.  If  Ura  is  approximated  by  and 
Equation  (B7)  is  substituted  into  Equation  (B6),  the  final 
expression  relating  a  to  0  is  given  by 


(n,  )2 


.  0  .  (m  +  1 ) 


(B8) 


Various  authors  have  tabulated  values  of  n  variation 

0 

with  8  .  In  this  work  values  from  Reference  18  were  used 
to  convert  R«j,crit  results  from  References  19  and  7  from  a 
8  dependence  to  a  dependence  on  A  .  These  results  are  given 
in  Table  B. 
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000 

000 

667 

429 

429 

332 


0.250 

0.176 

0.111 

0.111 

0.053 

0.053 

0.026 

0.000 

0.000 

•0.024 

■0.024 

•0.048 

-0.048 

•0.065 

-0.090 

■0.090 


3.372 

3.372 

3.564 

3.769 

3.769 

3.899 

4.045 

4.206 

4.393 

4.393 

4.626 

4.626 

4.767 

4.924 

4.924 

5.244 

5.244 

5.356 

5.356 

5.624 

6.771 

6.771 


11.370 

11.370 

8.468 

6.088 

6.088 

5.067 

4.091 

3.122 

2.144 

2.144 

1.126 

1.126 

0.583 

0.000 

0.000 

-0.671 

-0.671 

-1.366 

-1.366 

-2.069 

-4.162 

-4.162 


Rfi1 >crit 


12400 

12490 

10920 

8640 

8890 

7680 

6230 

4550 

2955 

2830 

1658 

1380 

865 

680 

520 

354 

318 

126 

199 

138 

0 

67 


S~J 


Boundary-Layer  Thicknesses 


For  the  two  internal  flow  cases  the  boundary-layer 
displacement  thicknesses  and  the  boundary-layer  moment 
thicknesses  were  approximated  following  a  method  given  by 
Dr.  J.  Hitchcock  (Ref.  20). 


Eacallfil  Elai  Elates 


The  common  definition  for  the  boundary-layer 
displacement  thickness  was  given  by: 


«1  = 


/ 


u(y ) 


(Cl) 


Here,  a  flow  area  was  equated  to  the  integration  of 

(1  -  u(y)/U)  over  a  differential  perimeter  area  in  the 


form : 


(flow  area)  = 


/ 


u(y) 


)  d(perimeter  area)  (C2) 


flow  area 


where  from  Figure  Cl  a 

(flow  area)  s  2  *•  w  •  51 
d  (perimeter  area)  s  2wdy  +  2hdz 


(C3) 

(C4) 


For  two  infinite  parallel  plates  (Fig.  Clb)  the  differential 
perimeter  area  became  2wdy.  Using  a  parabolic  distribution 


for  the  velocity, 
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The  expression  for  the  displacement  thickness  was  then 


h  2v 

*1  ■  /  <’  -  -H  +  ZE>  dy 

o 

or, 

«1  =  h/3 

2 .  Momentum  Xhisknsas 

The  momentum  thickness  is  commonly  given  as: 

52  =  f  2<I>  (1  .  )dy 

J  u  u 

o 

The  approximation  used  here  was 


or, 


(flow  area)  =  /  (2igi  C,  -  )) 

flow  area 

d  (perimeter  area) 


2w,2  =  /  *»>  (1  -  )  (2wdy ) 

o 


Again,  using  a  parabolic  velocity  distribution  yielded 


c 


B.  £1B£  flflH 


Thickness 


The  approximation  used  here  for  the  boundary-layer 
thickness  for  the  pipe  was 


r  U(r) 

(flow  area)  =  /  (1  -  - )  d(perimeter  area)  (C12) 

J  U 

o 


where  (see  Fig.  20 

flow  area  =  2itR6^  _  ir&|2 

perimeter  =  2irrdr 
area 

With  a  parabolic  velocity  distribution  of, 

U(r)  ,  r^ 

U  R2 

a  displacement  thickness  expression  was  found  to  be 


( C 1 3 ) 


.  r  p2 

2  R«i  -  6-|2  =  2R  J  (1  -  (1  -  — ))  rdr 


(Cl  4) 


2R6<j  -  5^  =  R/2 


(C15) 


This  expression  resulted  in  a  value  of  1  of  approximately 
(0.293)R. 


2 .  Mfifflfijitiiffl  Thickness 


The  boundary-layer  momentum  thickness  was  approximated 
using: 

2  R 

(2irR62  -  H6  (1  -  Hill)  (2*rdr)  (C16) 

J  V  U 

o 

Substituting  in  Equation  (C13)  yielded 

R 

2  2 

(2Rs ?  -  S22>  =  [  (1  «  ~)  (~)  4rdr  (C17) 

i  R2  R2 


or, 


2R62  -  «22 


r£ 

3 


which  when  solved  for  a  value  of  s2  =  (0.1835)R. 


V?l 

■,A\ 


m 


AbeshaI-x  u 
Test  Trials 


The  following  computer  program,  written  in  its  present 
form  by  J.  Lawrence  (Ref.  17),  was  used  to  test  the 
prediction  methods  developed  in  this  work.  Table  D  shows 
the  results  of  the  runs  made.  The  original  function  of  this 
program  was  to  investigate  the  effect  of  a  pitching  angle  of 
attack  on  the  flow  parameters,  Cl,  Cd,  etc.  Here  a  was  set 
at  zero. 


TABLE  D 


RL(X106) 


TRIAL 


CL 

n  i  =  const 

(X/l) 

It  S  I  \  # 

(X/l) 

1.7 

0.6 

0.190 

0.130 

1.7 

0.5 

0.195 

0.140 

1.7 

0.4 

0.209 

0.154 

1.7 

0.2 

0.229 

0.170 

1.7 

0.0 

0.245 

0.200 

1.7 

-0.2 

0.270 

0.217 

5.0 

0.4 

0.100 

0.124 

5.0 

0.2 

0.10$ 

0.135 

5.0 

0.0 

0.124 

0.174 

5.0 

-0.2 

0.132 

0.184 

aA.a^T  «J5  ■XmJlYV «_s  AA  A  AA 


PROGRAM  POHL3 

COMPLEX  CMPLX,EI,Z,DZETA,W 
OPEN  < 15, FILE-  FLOWIN' > 

REWIND  15 

OPEN  ( 16 ,FILE='docout ' ) 

EI=(0. , 1 ♦ ) 

RADIUS-1 *131 
AMU=- . 131 

READ  <15,*,END=i00)  ALPHA, ADOT1 ,UINF 

ALPH1-ALPHA 

CON-3.1415927/180. 

THETA-180. 

TIME-0.0 

CALL  DS< 180* , RADIUS, CON, AMU, XLE,YLE> 

CALL  DS ( 0 . 0 , RADIUS , CON , AMU , XTE , YTE > 

XLE-ABS(XLE) 

CHORD-XLE+XTE 

PITCH=AD0T1*C0N*0 . 5*CH0RD/UINF 

K-100 

Kl-K+1 

K2=(ALPHA+150)*100+K 

WRITE( 16,30 JUINF 
WRITE ( 16,40) ALPHA 
WRITE (16, 50 )AD0T1 
WRITE(16,52) 

ANGLE=ALPHA+THETA 

CALL  U ( ANGLE , RADIUS , CON , El , U I NF , AMU , ALPHA , UO ) 

CALL  DS ( ANGLE , RADIUS , CON , AMU , XO , YO  > 

ANGLE-ANGLE-0.01 

CALL  U < ANGLE , RADIUS , CON , E I , U I NF , AMU , ALPHA , U1 > 

CALL  DS ( ANGLE , R AD IUS , CON , AMU , X 1 , Y 1 > 

ANGLE- ANGLE-0. 01 

CALL  U < ANGLE , RADIUS , CON , E I , U INF , AMU , ALPHA , U2 ) 

CALL  DS( ANGLE, RADIUS, CON, AMU, X2,Y2) 

DS2*(SQRT ( (X2-X1 )##2+<Y2-Yl >##2) ) /CHORD 
DS1- ( SORT << Xl-XO ) **2+ < Y1 -YO ) **2 ) ) /CHORD 

Stagnation  point  velocity  gradient  computed  using  a 
forward  difference  method?  all  other  velocity  gradients 
computed  using  central  difference  method. 

DUDS- ( U2-U0 ) / <  DS 1 +DS2  > 

Second  derivative  of  velocity  computed  using  a 
Taylor's  Series  expansion. 

D2UDS2-(U2-2.#Ul+U0>/< (DSl+DS2>/2» >##2 


*********  *  *  *  * * 


Enter  initial  boundary  layer  parameters. 

RLAMDA=7 . 052 
RK=0 .0770 
FK=0 .0 

DZDS=-0 . 0652*D2UDS2/  <  DUDS**2  ) 

ZZ=RK/DUDS 

N=50 

ANGLE® ALPHA+THETA-0 .01 
XOC®  <  XO+XLE )  /CHORD 
F2K=0  ♦  C 

WRITE (16,1) X0C , UO , F2K , RLAMDA , FK , RK , ZZ , DZDS 

AD0T=0.0 
DO  10  J®1 ,  K 

Function  of  this  loop  is  to  compute  boundary  layer 
parameters  at  stagnation  point*  allowing  the 
boundary  layer  to  steady-out  before  subjecting  it 
to  a  pitching  airflow. 

N=N+1 

Compute  pertinent  boundary  layer  porometers. 

ZZ=DZDS*DS1+ZZ 

RK»ZZ*DUDS 

FN=.47-6.*RK 

DZDS®FK/U1 

DELT=CH0RD*DS1/U1 

TIME=TIME+DELT 

CALL  U  <  ANGLE , RAD I US , CON  , E I , UINF , AMU  *  ALPHA , U2 ) 

DUDT® (U2-U1 )/DELT 
ANGLE®ANGLE-0 .01 
ANGLEl®ANGLE-0 .01 

CALL  U  <  ANGLE1 *  RADIUS , CON , El  *  UINF  * AMU  *  ALPHA  *  U2 ) 

CALL  DS( ANGLE1 , RADIUS, CON, AMU, X2,Y2> 

ANGLE0=ANGLE+0 .01 

CALL  U  <  ANGLEO  *  RADIUS , CON  *  El , UINF  , AMU  *  ALPHA , UO ) 

CALL  DS( ANGLEO, RADIUS  ,C0N, AMU,X0, YO) 

CALL  U <  ANGLE , RAD I US , CON , E I , U INF , AMU , ALPHA , U1 > 

CALL  DS < ANGLE , RADIUS , CON , AMU , XI , Y1 ) 

DS1=<SGRT<  <Xl-X0)**2+< Yl-Y0>**2> >/CH0RD 
DS2® ( SORT ( <  X2-X1 ) #*2  +  <  Y2-Y1 ) **2 ) ) /CHORD 


DSS»DS1+DS2 
DUDS® ( U2-U0 ) /DSS 
UDUDS»U1*DUDS 
XOC® <  XI +XLE) /CHORD 


***  *****  ***  ***** 


I~<N.L7.50>  GO  70  10 
N-0 

WRITE <16, 1)X0C,U1,F2K,RLAMDA,FK,RK,ZZ,DZDS 

10  CONTINUE 

* 

AD0T=AD0T1 
N*0 

00  20  J=K1 »N2 

Function  of  this  loop  is  to  compute  the  behavior 
of  the  boundary  layer  as  it  is  subjected  to  a 
pitching  airfoil* 

N-N+l 

Compute  the  pertinent  boundary  layer  parameters* 

ZZ»DZDS*DS1+ZZ 
RK*ZZ*  <  DUDS+DUDT/U1 ) 

CALL  POHL<RK,RLAMDA) 

DEL2=37 . /31 5 . -RLAMDA/945  *  - <  RLAMDA##2 ) /9072 ♦  ■ 

FK*2  *  *DEL2* ( 2  *  - . 3683*RLAMDA+ . 0104*RLAMDA**2+ 

+  <RLAMDA**3)/4536> 

F2K*< . 3-RLAMDA/120* )/DEL2 
DZDS= <  FK+ <  4 . +F2K ) *ZZ*DUDT/U1 > /U1 

Compute  the  time  increment  for  a  particle  to 
travel  from  point  <i>  to  point  <i+l>» 

DELT«CH0RD*DS1/U1 
TIME-TIME+DELT 
DALPHA*DELT*AD0T 
ANGLE»ANGLE+DALPHA 
ALPH1 *ALPH1 +DALPHA 

CALL  U < ANGLE , RADIUS , CON , E I , UINF , AMU , ALPH 1 , U2 ) 

Compute  the  unsteady  velocity  gradient* 

DUDT*<U2-U1 )/DELT 
ANGLE*ANGLE-0 .01 
ANGLE1* ANGLE-0. 01 

CALL  U < ANGLE1 , RADIUS , CON , El , UINF , AMU , ALPH1 , U2 ) 
CALL  DS<ANGLE1, RADIUS, CON, AMU, X2,Y2> 
ANGLE0«ANGLE+0 .01 

CALL  U < ANGLEO , RADIUS , CON , E I , U I NF , AMU , ALPH1 , UO ) 
CALL  DS < ANGLEO , RADIUS , CON , AMU ,X0 , YO ) 

CALL  U< ANGLE, RADIUS, CON, El, UINF, AMU, ALPH1,U1) 

‘  CALL  DS<ANGLE, RADIUS, CON, AMU, XI, Yl> 

* 

*  Compute  arc  length  and  velocity  gradient* 

* 
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CS 1  =  (  SORT  «  i  X  J  -/ C  U  *2  +  ( Y 1  -YO )  **2  >'  )'7CH0RD 

DS2- < SQRT  £ (X2-Xl)»*2r(Y2-Yl)**2) > /CHORD 

DSS*DS1+DS2 

DUDS- ( U2-U0  > /DSS 

UDUDS=U1*DUDS 

XOC=< XI +XLE) /CHORD 

Stop  the  computation  at  the  quarter-chord* 

IF (XOC«GE *0*500)  GO  TO  25 
IF<N.LT*250)  GO  TO  20 
N*0 

WRI TE  < 1 6 , 1 ) XOC ,  U1 , F2K , RLAMDA , FK , RK , ZZ , DZDS 
20  CONTINUE 

25  WRITE (16*1) XOC  *  U1 » F2K » RLAMDA *  FK  * RK*ZZ  * DZDS 
WRITE ( 16*45) ALPH1 
WRITE < 16 *55) PITCH 
WRITE<16,60)RK 
WRITE( 16*80) 

WRITE(16,S1)TIME 

1  FORMAT  (4X»F6«3*2(4X»F10*3)*4X»F7*3*4X»F7.4*4X*F8*4*2(4X»E9. 

30  FORMAT < 1H1 * 'BOUNDARY-LAYER  PARAMETERS  FOR  * ,F6.2* *FT/SEC*/) 

40  FORMAT  <  *  INITIAL  ANGLE  OF  ATTACK:  *,F6*3**  DEGREES*/) 

45  FORMAT </*  FINAL  ANGLE  OF  ATTACK:  **F6,3**  DEGREES*/) 

50  FORMAT <  *  PITCH  RATE:  *,F7,3,*  DEGREES/SEC*/) 

52  FORMAT (6X* *XOC* * 10X , *U* * 1 IX * ’ F2K * * 9X * 'LAMBDA* *8X* *FK* * 

+  9X**RK,*11X»*Z,»11X**  DZDS  * / ) 

55  FORMAT (*  PITCH  PARAMETER:  *,F7.5/) 

60  FORMAT C  K  AT  THE  QUARTER-CHORD:  **F8*4/) 

80  FORMAT <*  TIME  TO  REACH  THE  QUARTER-*) 

81  FORMAT  <  *  CHORD  FROM  THE  STAGNATION  POINT:  *,F7*5,*  SEC*/) 

STQP 

00  END 


SUBROUTINE  U < ANGLE * RADIUS * CON * El * UINF * AMU , ALPHA * UU ) 
COMPLEX  CMPLX*Z*EI*DZETA,W 

Function  of  this  subroutine  is  to  to  compute  the  local 
value  of  velocity  on  a  JouKowsKi  Airfoil  using  complex 
potential  flow  theory* 

X»RAD1US*C0S < ANGLE*CON ) 

Y*RADIUS#SIN ( ANGLE*CON> 

Z»CMPLX(X,Y) 

W*UINF*< (1**0. )-<RADIUS**2)/Z**2+ 
f  <2.*EI*RADIUS*SIN(ALPHA*C0N) )/Z) 

X=X+AMU 


*•  *  •*  * 


Z  changed  to  repress;.  L  vula-s  coordinates  used  in 

the  transformation  equation* 

Z=CMPLX(X»Y) 

DZETA= <  Z**2- ( RAOIUS+AMU ) **2 ) /Z**2 
UU*CABS<W>/CABS<DZETA> 

RETURN 
END 

* 

« 

SUBROUT I NE  DS < ANGLE , RADI US , CON , AMU ,X»Y> 

.  COMPLEX  CMPLXt Z 

* 

*  Function  of  this  subroutine  is  to  compute  the  arc 

*  length  between  two  points  on  the  JoukowsKi  Airfoil* 

* 

X=RADIUS*COS  <  ANGLE*CON ) +AMU 
Y=RADIUS*SIN ( ANGLE*CON ) 

Z*CMPLX<X,Y> 

Z=Z+<  <RADIUS+AMU)**2)/Z 
X-REAL(Z) 

Y*AIMAG(Z> 

RETURN 

END 

* 

* 

SUBROUTINE  POHL  <  RK  t RLAMDA ) 

* 

*  Function  of  this  subroutine  is  to  compute  the  value 

*  of  the  separation  parameter T  Lamda r  given  a  value 

*  of  Kf  as  computed  in  the  main  program. 

* 

RKJ  =- . 160 

RK2*-.112 

RK3=0.00 

RK4=0 .06 

RK5*0.076 

RK6=0.086 

RK7*0.0949 

IF(RK.LE.RKl)  GO  TO  10 
IF<RK.LE*RK2)  GO  TO  20 
IFtRK.LE  .RK3)  GO  TO  30 
IF (RK.LE *RK4)  GO  TO  40 
IF<RK.LE.RK5>  GO  TO  50 
IF (RK.LE.RK6)  GO  TO  60 
IF(RK»GT *RN7)  GO  TO  70 

* 

RLAMDA* .0149**2- ( RN-0 . 08 ) **2 
RL AMDA* 12.-1 00 »*SGRT (RLAMDA) 

RETURN 


RLAMDA—  ( 2 .  /  ♦01^..'  ♦r-.i\  r  i  .  0 

RETURN 

RLAMTiA=  <  4 »  / .  044 )  #RK+2  *  18 
RETURN 

RLAMDA=<10./.14)*RK 

RETURN 

RLAMDA=83.33*RK 

RETURN 

RLAMDA=-1 .9+115. *RK 
RETURN 

RLAMDA=~6 . 54+176 . *RK 
RETURN 

RL.AMDA=12 . 

RETURN 
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